/********************************************************
* filename : tms2dwt.c					*
* startdate : Thu, November 5 2009			*
* lastupdate : Thu, October 14 2010			*
* author : Limthong K. Kriangkrai			*
* location : NII Tokyo, Japan				*
*							*
* Convert a time series file to a wavelet file		*
********************************************************/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_wavelet.h>
#include <unistd.h>
#include <ctype.h>

#define MAXLINE 300 //Maximum line of data in the time series file.

int main(int argc, char *argv[]) {
	//To check arguments.
	int dflag = 0;
	int level = 1;
	char *ofile = NULL;
	int n = 131072;
	char dwtAfile[100];
	char dwtDfile[100];

	int c;
	opterr = 0;

	if(argc < 2) {
		printf("Usage: %s TMSFILE [OPTION]\n", argv[0]);
		printf("Convert a time series file to a wavelet file\n\n");
		printf("  -d Show the detail part.\n");
		printf("  -j NUM convert to wavelet level NUM (default 1).\n");
		printf("  -o FILE write output to FILE.\n");
		printf("  -s NUM set number of wavelet data to NUM (default 131072).\n\n");
		printf("Report bugs to <krngkr@nii.ac.jp>.\n");
		return 1;
	}
	char *tmsfile = argv[1];
	//To check file existing.
	FILE *testFile;
	testFile = fopen(tmsfile, "r");
	if(testFile == NULL) {
		perror ("tms2dwt");
		return 1;
	}
	while((c = getopt(argc, argv, "dj:o:s:")) != -1)
		switch(c) {
			case 'd':
				dflag = 1;
				break;
			case 'j':
				level = atoi(optarg);
				break;
			case 'o':
				ofile = optarg;
				break;
			case 's':
				n = atoi(optarg);
				break;
			case '?':
				if(optopt == 'o' || optopt == 'j' || optopt == 's')
					fprintf(stderr, "Option -%c requires an argument.\n", optopt);
				else if(isprint(optopt))
					fprintf(stderr, "Unknown option `-%c'.\n", optopt);
				else
					fprintf(stderr, "Unknown option character `\\x%x'.\n", optopt);
				return 1;
		}
	//To check DWT level.
	if(level <= 0) {
		printf("DWT Level %d : Can't decomposition.\n", level);
		return 1;
	} else if(level > (log(n)/log(2))) {
		printf("Maximum of DWT level is %g\n", log(n)/log(2));
		return 1;
	}
	//End check arguments.

	if(ofile != NULL) {
		strcpy(dwtAfile,ofile); strcat(dwtAfile,".A");
		strcpy(dwtDfile,ofile); strcat(dwtDfile,".D");
		printf("%s > %s, %s\n", tmsfile, dwtAfile, dwtDfile);
	}

	char *packettime;
	int timestamp;
	double packet, size, flow, saddr, daddr, sport, dport, dtaddr, dtport;

	char  packettimearray[n][20];
	double *data = malloc(n * sizeof(double));
	double *datacA = malloc(n * sizeof(double));
	double *datacD = malloc(n * sizeof(double));
	double *dataA = malloc(n * sizeof(double));
	double *dataD = malloc(n * sizeof(double));

	double *apacket = malloc(n * sizeof(double));
	double *asize = malloc(n * sizeof(double));
	double *aflow = malloc(n * sizeof(double));
	double *asaddr = malloc(n * sizeof(double));
	double *adaddr = malloc(n * sizeof(double));
	double *asport = malloc(n * sizeof(double));
	double *adport = malloc(n * sizeof(double));
	double *adtaddr = malloc(n * sizeof(double));
	double *adtport = malloc(n * sizeof(double));

	double *dpacket = malloc(n * sizeof(double));
	double *dsize = malloc(n * sizeof(double));
	double *dflow = malloc(n * sizeof(double));
	double *dsaddr = malloc(n * sizeof(double));
	double *ddaddr = malloc(n * sizeof(double));
	double *dsport = malloc(n * sizeof(double));
	double *ddport = malloc(n * sizeof(double));
	double *ddtaddr = malloc(n * sizeof(double));
	double *ddtport = malloc(n * sizeof(double));

	//Initial data to zero.
	for(timestamp = 0; timestamp < n; timestamp++) {
		strcpy(packettimearray[timestamp],"NULL");
		data[timestamp] = 0;
		apacket[timestamp] = 0; dpacket[timestamp] = 0;
		asize[timestamp] = 0; dsize[timestamp] = 0;
		aflow[timestamp] = 0; dflow[timestamp] = 0;
		asaddr[timestamp] = 0; dsaddr[timestamp] = 0;
		adaddr[timestamp] = 0; ddaddr[timestamp] = 0;
		asport[timestamp] = 0; dsport[timestamp] = 0;
		adport[timestamp] = 0; ddport[timestamp] = 0;
		adtaddr[timestamp] = 0; ddtaddr[timestamp] = 0;
		adtport[timestamp] = 0; ddtport[timestamp] = 0;
	}

	int columnindex;

	//Read a time series file.
	char line[MAXLINE];
	char *sub_string;
	FILE *tFile;
	if((tFile = fopen(tmsfile, "r")) != NULL) {//Open the time series file for read data.
		while(fgets(line, MAXLINE, tFile) != NULL) {
			packettime = strtok(line, ",");
			//printf("%s, ", packettime);
			columnindex = 0;
			while((sub_string = strtok(NULL, " \t\n")) != NULL) {
				switch(columnindex) {
					case 0:
						timestamp = atoi(sub_string);
						//printf("%d\t", timestamp);
						break;
					case 1:
						packet = atof(sub_string);
						//printf("%.10g ", packet);
						break;
					case 2:
						size = atof(sub_string);
						//printf("%.10g ", size);
						break;
					case 3:
						flow = atof(sub_string);
						//printf("%.10g ", flow);
						break;
					case 4:
						saddr = atof(sub_string);
						//printf("%.10g ", saddr);
						break;
					case 5:
						daddr = atof(sub_string);
						//printf("%.10g ", daddr);
						break;
					case 6:
						sport = atof(sub_string);
						//printf("%.10g ", sport);
						break;
					case 7:
						dport = atof(sub_string);
						//printf("%.10g ", dport);
						break;
					case 8:
						dtaddr = atof(sub_string);
						//printf("%.10g ", dtaddr);
						break;
					case 9:
						dtport = atof(sub_string);
						//printf("%.10g\n", dtport);
						break;
				}
			columnindex++;
			}
			strcpy(packettimearray[timestamp],packettime);
			apacket[timestamp] = packet;
			asize[timestamp] = size;
			aflow[timestamp] = flow;
			asaddr[timestamp] = saddr;
			adaddr[timestamp] = daddr;
			asport[timestamp] = sport;
			adport[timestamp] = dport;
			adtaddr[timestamp] = dtaddr;
			adtport[timestamp] = dtport;
		}
	}
	fclose(tFile);//Close a time series file.

	//Initial wavelet variable.
	int nd;
	gsl_wavelet *w;
	gsl_wavelet_workspace *work;
	//w = gsl_wavelet_alloc(gsl_wavelet_daubechies, 4);//Can change the mother of wavelet here.
	w = gsl_wavelet_alloc(gsl_wavelet_haar, 2);
	work = gsl_wavelet_workspace_alloc(n);

	int i, j, k;

	for(k = 0; k < 9; k++){
		nd = n;//Reset nd to n.
		switch(k) {
			case 0:
				for(i = 0; i < n; i++) {
					data[i] = apacket[i];
				}
				break;
			case 1:
				for(i = 0; i < n; i++) {
					data[i] = asize[i];
				}
				break;
			case 2:
				for(i = 0; i < n; i++) {
					data[i] = aflow[i];
				}
				break;
			case 3:
				for(i = 0; i < n; i++) {
					data[i] = asaddr[i];
				}
				break;
			case 4:
				for(i = 0; i < n; i++) {
					data[i] = adaddr[i];
				}
				break;
			case 5:
				for(i = 0; i < n; i++) {
					data[i] = asport[i];
				}
				break;
			case 6:
				for(i = 0; i < n; i++) {
					data[i] = adport[i];
				}
				break;
			case 7:
				for(i = 0; i < n; i++) {
					data[i] = adtaddr[i];
				}
				break;
			case 8:
				for(i = 0; i < n; i++) {
					data[i] = adtport[i];
				}
				break;
		}

		//Start loop level.
		for(i = 0; i < level; i++) {
			gsl_wavelet_transform_forward (w, data, 1, nd, work);

			for(j = 0; j < nd; j++) {
				if(j < nd/2) {
					datacA[j] = data[j];
				} else {
					datacA[j] = 0;
				}
				if(j >= nd/2) {
					datacD[j] = data[j];
				} else {
					datacD[j] = 0;
				}
			}

			//gsl_wavelet_transform_inverse(w, count, 1, nd, work);
			gsl_wavelet_transform_inverse(w, datacA, 1, nd, work);
			gsl_wavelet_transform_inverse(w, datacD, 1, nd, work);

			for(j = 0; j < nd; j++) {
				if(j < nd/2) {
					data[j] = datacA[2*j];
				} else {
					data[j] = 0;
				}
			}
			nd = nd/2;

		}//End loop level.

		int step = pow(2,level-1);
		for(i = 0; i < n/step; i++) {
			for(j = 0; j <= pow(2,level-1); j++) {
				dataA[(i*step)+j] = datacA[i];
				dataD[(i*step)+j] = datacD[i];
			}
		}

		switch(k) {
			case 0:
				for(i = 0; i < n; i++) {
					apacket[i] = dataA[i];
					dpacket[i] = dataD[i];
				}
				break;
			case 1:
				for(i = 0; i < n; i++) {
					asize[i] = dataA[i];
					dsize[i] = dataD[i];
				}
				break;
			case 2:
				for(i = 0; i < n; i++) {
					aflow[i] = dataA[i];
					dflow[i] = dataD[i];
				}
				break;
			case 3:
				for(i = 0; i < n; i++) {
					asaddr[i] = dataA[i];
					dsaddr[i] = dataD[i];
				}
				break;
			case 4:
				for(i = 0; i < n; i++) {
					adaddr[i] = dataA[i];
					ddaddr[i] = dataD[i];
				}
				break;
			case 5:
				for(i = 0; i < n; i++) {
					asport[i] = dataA[i];
					dsport[i] = dataD[i];
				}
				break;
			case 6:
				for(i = 0; i < n; i++) {
					adport[i] = dataA[i];
					ddport[i] = dataD[i];
				}
				break;
			case 7:
				for(i = 0; i < n; i++) {
					adtaddr[i] = dataA[i];
					ddtaddr[i] = dataD[i];
				}
				break;
			case 8:
				for(i = 0; i < n; i++) {
					adtport[i] = dataA[i];
					ddtport[i] = dataD[i];
				}
				break;
		}
	}

	FILE *dAFile, *dDFile;
	if(ofile != NULL) {
		dAFile = fopen(dwtAfile,"w");//Open the approximation wavelet file.
		dDFile = fopen(dwtDfile,"w");//Open the detail wavelet file.
	}
	for(i = 0; i < 86400; i++) {
		if(strcmp(packettimearray[i],"NULL") != 0) {
			if(ofile == NULL && dflag == 0) printf("%s, %d\t %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g\n", packettimearray[i], i, apacket[i], asize[i], aflow[i], asaddr[i], adaddr[i], asport[i], adport[i], adtaddr[i], adtport[i]);
			else if(ofile == NULL && dflag == 1) printf("%s, %d\t %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g\n", packettimearray[i], i, dpacket[i], dsize[i], dflow[i], dsaddr[i], ddaddr[i], dsport[i], ddport[i], ddtaddr[i], ddtport[i]);
			else if(ofile != NULL) {
				fprintf(dAFile, "%s, %d\t %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g\n", packettimearray[i], i, apacket[i], asize[i], aflow[i], asaddr[i], adaddr[i], asport[i], adport[i], adtaddr[i], adtport[i]);
				fprintf(dDFile, "%s, %d\t %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g\n", packettimearray[i], i, dpacket[i], dsize[i], dflow[i], dsaddr[i], ddaddr[i], dsport[i], ddport[i], ddtaddr[i], ddtport[i]);
			}
			//else if(dwtfile != NULL && dflag == 0) fprintf(dFile, "%s, %d\t %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g\n", packettimearray[i], i, apacket[i], asize[i], aflow[i], asaddr[i], adaddr[i], asport[i], adport[i], adtaddr[i], adtport[i]);
			//else if(dwtfile != NULL && dflag == 1) fprintf(dFile, "%s, %d\t %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g\n", packettimearray[i], i, dpacket[i], dsize[i], dflow[i], dsaddr[i], ddaddr[i], dsport[i], ddport[i], ddtaddr[i], ddtport[i]);
		}
	}
	if(ofile != NULL) {
		fclose(dAFile);//Close the approximation wavelet file.
		fclose(dDFile);//Close the detail wavelet file.
	}

	gsl_wavelet_workspace_free(work);
	gsl_wavelet_free(w);

	free(data);
	free(datacA);
	free(datacD);
	free(dataA);
	free(dataD);

	return 0;
}
